package org.chorus;

import jwo.landserf.gui.GISFrame;
import jwo.landserf.gui.SimpleGISFrame;
import jwo.landserf.process.io.FileHandler;
import jwo.landserf.process.io.LandSerfIO;
import jwo.landserf.structure.ColourTable;
import jwo.landserf.structure.RasterMap;
import jwo.landserf.structure.SpatialModel;
import processing.core.PApplet;
import org.chorus.utils.LUDecomp;

//  **********************************************************************************
/** Sketch for detrending a raster surface. Taking a raster as input, the quadratic 
 *  polynomial that best fists the surface is found through least squares and the 
 *  residual values are saved as a new raster. This can be a useful pre-processing 
 *  step when performing variogram and PDE analysis.
 *  @author Jo Wood
 *  @version 1.2.1, 13th March, 2013
 */
//  **********************************************************************************

/* This file is part of the Chorus river bed analysis software package. Chorus is free software:
 * you can redistribute it and/or modify it under the terms of the GNU General Public License
 * as published by the Free Software Foundation, either version 3 of the License, or (at your
 * option) any later version.
 * 
 * Chorus is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
 * without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  
 * See the GNU General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License along with this
 * source code (see COPYING included with this source code). If not, see 
 * http://www.gnu.org/licenses/.
 */

@SuppressWarnings("serial")
public class Detrender extends PApplet
{
	// ------------------------ Object and class variables --------------------------
	
	private String statusMessage;				// Used for error reporting to screen.
	private static String propertiesFile;		// Name and location of Chorus properties file.
	private static final String VERSION = "1.2.1, 13th March, 2013";
	
	// ------------------------------ Starter method -------------------------------

	/** Starts the sketch as an application.
	 *  @param args Command line arguments (ignored). 
	 */
	public static void main(String[] args)
	{   
		if (args.length > 0)
		{
			propertiesFile = args[0];
		}
		else
		{
			propertiesFile = "scripts/chorus.props";
		}
		PApplet.main(new String[] {"org.chorus.Detrender"});
	}

	// ---------------------------------- Methods ----------------------------------
	

	/** Initialises the detrender sketch.
	 */
	public void setup()
	{
		size(800,200);
		smooth();
		textFont(createFont("GeosansLight.otf",18));
		
		statusMessage = "Problem detrending files";

		if (frame !=null)
		{
			frame.setTitle("Surface Detrender "+VERSION);
		}
		
		// Get the list of files to detrend from the properties file.
		PropertyHandler props = new PropertyHandler(propertiesFile,this);
		
		if (props.foundPropertiesFile() == false)
		{
			statusMessage = "No properties file found";
		}
		else
		{
			String baseDir = props.getBaseDir();
			String[] rastNames = props.getDetrendNames();
			
			if (rastNames.length == 0)
			{
				statusMessage = "No files to detrend found in "+propertiesFile;
			}
			else
			{
				GISFrame gisFrame = new SimpleGISFrame();
				
				for (String rastName : rastNames)
				{
					String rastBaseName = rastName;
					String extension = ".srf";
					int dotLocation = rastName.lastIndexOf(".");
					if (dotLocation >0)
					{
						// File extension provided.
						rastBaseName = rastBaseName.substring(0,dotLocation);
						extension = rastName.substring(dotLocation);
					}

					int rastType = FileHandler.guessFileType(rastBaseName+extension);
					RasterMap origRaster = (RasterMap)FileHandler.readFile(baseDir+rastBaseName+extension, rastType, gisFrame, null);
										
					if (origRaster == null)
					{
						statusMessage = "Cannot read raster from "+baseDir+rastBaseName+extension;
						return;
					}
					detrend(origRaster,baseDir+rastBaseName+"Detrended.srf");
				}
				// If we get to this point, we must have read and detrended all files successfully.
				statusMessage = (rastNames.length>1?"DEMs":"DEM")+" detrended";
			}
		}
	}

	/** Draws the user interface for the detrender.
	 */
	public void draw()
	{
		background(255);

		fill(0);
		textSize(22);
		textAlign(CENTER,CENTER);
		text(statusMessage,10,10,width-10,height-10);
		noLoop();
	}

	// ------------------------------ Private methods ------------------------------

	/** Detrends the given raster and saves the detrended version with the given file name.
	 *  @param origRaster Raster to detrend.
	 *  @param rasterFileName File name of new raster that will contain the detrended surface.
	 */
	private void detrend(RasterMap origRaster, String rasterFileName)
	{
		double cx = origRaster.getNumCols()/2f;
		double cy = origRaster.getNumRows()/2f;

		double[][] normal = new double[6][6];
		findNormal(origRaster,cx,cy,normal);
		LUDecomp luDecomp = new LUDecomp(normal);    // Decompose matrix.

		// Find the observed values and place in vector.
		double obs[] = new double[6]; // Observed value vector.

		for (int row=0; row<origRaster.getNumRows(); row++)
		{
			for (int col=0; col<origRaster.getNumCols(); col++)
			{   
				double z = origRaster.getAttribute(row, col);
				if (z != SpatialModel.NO_VALUE)
				{
					// Create local coordinates.
					double x = col-cx;
					double y = row-cy;

					obs[0] += z*x*x;
					obs[1] += z*y*y;
					obs[2] += z*x*y;
					obs[3] += z*x;
					obs[4] += z*y;
					obs[5] += z;
				}
			}
		}

		// Solve equations based on observed vector.
		luDecomp.luBackSub(obs);

		// Create new surface based on difference between original surface and quadratic trend.
		RasterMap newRaster = new RasterMap(origRaster);
		for (int row=0; row<newRaster.getNumRows(); row++)
		{
			for (int col=0; col<newRaster.getNumCols(); col++)
			{       
				double origZ = origRaster.getAttribute(row, col);
				if (origZ != SpatialModel.NO_VALUE)
				{
					double x = col-cx;
					double y = row-cy;
					double trend = obs[0]*x*x + obs[1]*y*y + obs[2]*x*y + obs[3]*x +obs[4]*y + obs[5];
					newRaster.setAttribute(row, col, (float)(origZ-trend));
				}
			}
		}

		// Save new raster.
		newRaster.getHeader().setTitle("Detrended "+origRaster.getHeader().getTitle());
		newRaster.getHeader().setNotes("Detrended using quadratic function z = "+
				                       obs[0]+"x^2 + "+obs[1]+"y^2 + "+obs[2]+"xy + "+obs[3]+"x + "+obs[4]+"y + "+obs[5]+"\n"+
				                       origRaster.getHeader().getNotes());
		newRaster.setColourTable(ColourTable.getPresetColourTable(ColourTable.IMHOF_L3,newRaster.getMinAttribute(),newRaster.getMaxAttribute()));
		LandSerfIO.write(newRaster,rasterFileName);
	}


	/** Calculates normal equations based on a given window size.
	 */
	private void findNormal(RasterMap raster, double cx, double cy, double[][] normal)
	{
		double  x,y,     // Cross product parameters.
		x1  =0F, y1 =0F,
		x2  =0F, y2 =0F,
		x3  =0F, y3 =0F,
		x4  =0F, y4 =0F,
		xy2 =0F, x2y=0F,
		xy3 =0F, x3y=0F,
		x2y2=0F, xy =0F,
		N=0F;

		// Create the sums-of-squares and cross products matrix 
		for (int row=0; row <raster.getNumRows(); row++)
		{
			for (int col=0; col<raster.getNumCols(); col++)
			{
				if (raster.getAttribute(row, col) != SpatialModel.NO_VALUE)
				{
					x = col-cx;    // Convert to local grid 
					y = row-cy;    // coordinates.

					x4   += x*x*x*x;
					x2y2 += x*x*y*y;
					x3y  += x*x*x*y;
					x3   += x*x*x;
					x2y  += x*x*y;
					x2   += x*x;

					y4   += y*y*y*y;
					xy3  += x*y*y*y;
					xy2  += x*y*y;
					y3   += y*y*y;
					y2   += y*y;

					xy   += x*y;

					x1   += x;

					y1   += y;

					N    +=  1;
				}
			}
		}

		// Store cross-product matrix elements.
		normal[0][0] = x4;
		normal[0][1] = normal[1][0] = x2y2;
		normal[0][2] = normal[2][0] = x3y;
		normal[0][3] = normal[3][0] = x3;
		normal[0][4] = normal[4][0] = x2y;
		normal[0][5] = normal[5][0] = x2;

		normal[1][1] = y4;
		normal[1][2] = normal[2][1] = xy3;
		normal[1][3] = normal[3][1] = xy2;
		normal[1][4] = normal[4][1] = y3;
		normal[1][5] = normal[5][1] = y2;

		normal[2][2] = x2y2;
		normal[2][3] = normal[3][2] = x2y;
		normal[2][4] = normal[4][2] = xy2;
		normal[2][5] = normal[5][2] = xy;

		normal[3][3] = x2;
		normal[3][4] = normal[4][3] = xy;
		normal[3][5] = normal[5][3] = x1;

		normal[4][4] = y2;
		normal[4][5] = normal[5][4] = y1;

		normal[5][5] = N;
	}
}